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I. INTRODUCTION 

Lately there has been a renewed interest in the study 
of magnetic field induced quantum phase transitions in 
spin-one magnets with strong single-ion and exchange 
anisotropics^ ^. The discovery of S = 1 compounds, 
such as Y2BaNi05 or the organo-metallic frameworks 
[Ni(C2H8N2)2(N02)]C104 (NENP), [Ni(C2H8N2)2Ni(CN)4] 
(NENC) and [NiCl2-4SC(NH2)2] (DTN), fuelled experimen- 
tal and theoretical studies of the role of dimensionality and 
single-ion anisotrop>^^^^^. In most of the known S = 1 mag- 
nets, the ubiquitous Heisenberg exchange is complemented 
by single-ion anisotropy. The interplay between these in- 
teractions with external magnetic field and lattice geome- 
try can result in a rich variety of quantum phases and phe- 
nomena, including the Haldane phase of quasi- ID systems^^, 
field induced Bose Einstein condensation (EEC) of magnetic 
state and field induced ferronematic orderingEl, Interest in 
5* = 1 Heisenberg antiferromagnets with uniaxial exchange 
and single-ion anisotropics has gained additional impetus re- 
cently after it was shown to exhibit the spin analog of the elu- 
sive supersolidphase on a lattice over a finite range of mag- 
netic fields.C^ 

In contrast to its classical counterpart (S oo), S = 1 
systems become quantum paramagnets (QPM) for sufficiently 
strong easy-plane single-ion anisotropy. In other words, they 
do not order down to zero temperature, T = 0, because the 
dominant anisotropy term, D ^^.{S^Y > 0)' forces each 
spin to be predominantly in the non-magnetic l^*^ = 0) state: 
{S^ = 0\S!;\S^ = 0) = for = {x, y, z}. The application 
of a magnetic field , H, along the 2;-axis reduces the spin gap 
linearly in H since the field couples to a conserved quantity 
(total magnetization along the z-axis). The gap is closed at a 
quantum critical point (QCP) where the bottom of the = 1 
branch of magnetic excitations touches zero. This QCP be- 
longs to the EEC universality class and the gapless mode of 
low-energy = 1 excitations remains quadratic for small 
momenta, cx /c^ , because the Zeeman term commutes with 
the rest of the Hamiltonian. Since the dynamical exponent is 



z = 2, the effective dimension is d + 2 and the upper critical 
dimension is = 2. This, and analogous field-driven transi- 
tions, have been widely studied experimentally to demonstrate 
EEC related phenomena in many quantum magnets .'^J^^^^^t^U 
One of these magnets is the metal-organic framework DTN 
that we mentioned abovJ^Hl. 

The starting point of any theoretical study of a magnetic 
field induced phase transition in a QPM is to determine the 
Hamiltonian parameters, i.e., the exchange constants and the 
amplitude of the different anisotropics. The simplest way of 
extracting these parameters is to fit the branches of magnetic 
excitations that are measured with inelastic neutron scatter- 
ing (INS). The reliability of this procedure is normally lim- 
ited by the accuracy of the approach that is used to com- 
pute the dispersion relation of magnetic excitations. Numer- 
ical methods like Quantum Monte Carlo (QMC) and Den- 
sity Matrix Renormalization Group (DMRG) are very accu- 
rate, but they can only be applied under special circumstances. 
While the DMRG method^"^ has evolved to the extent that dy- 
namical properties such as the frequency and momentum de- 
pendence of the magnetic structure factor can be computed 
very accuratel>^^ , its application is restricted to quasi-one- 
dimensional magnets such as HPIP-CuEr4^^. On the other 
hand, QMC methods can only be applied to systems that have 
no frustration in the exchange interaction, i.e., that are free of 
the infamous sign problem. Consequently, it is necessary to 
find simple analytical approaches that are accurate enough to 
quantitatively reproduce the quantum phase diagram and the 
dispersion of magnetic excitations. 

One of the purposes of this work is to test different analyt- 
ical approaches against the results of accurate QMC simula- 
tions of a spin-one Heisenberg Hamiltonian with easy-plane 
single-ion anisotropy. The model is defined either on a square 
or on a cubic lattice to avoid frustration and make the QMC 
method applicable. Eesides being relevant for describing real 
quantum magnets, such as DTN, this model provides one of 
the simplest realizations of quantum paramagnetism and is 
ideal for testing methods that can be naturally extended to 
more complex systems. 
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II GENERALIZED SPIN WAVE APPROACH 



The generic S = 1 Heisenberg model with uniaxial single- 
ion anisotropy on an isotropic hyper-cubic lattice is given by 
the Hamiltonian: 

where the sum in the first term runs over nearest neighbor 
pairs {r^r'). D is the strength of the single ion-anisotropy, 
J is the exchange constant and hz = g/J^sH, where g is the 
g-factor and /i^ is the Bohr magneton . Henceforth, J is set to 
unity and all the parameters are expressed in units of J. In this 
work, we shall only consider models with spatially isotropic 
interactions, although the formalism can be straightforwardly 
generalized to anisotropic lattices. 

The {D^hz) quantum phase diagram of Hh is well known 
from mean field analysi^^^^, series expansion studied and 
numerical simulations ^"2. The D term splits the local spin 
states into = and = ±1 doublet. As we explained 
above, the ground state is a quantum paramagnet for large 
D 1, i.e., it has no long range magnetic order and there 
is a finite energy gap to spin excitations. At finite magnetic 
fields, the Zeeman term lowers the energy of the = 
state until the gap closes at a critical field he. A canted an- 
tiferromagnetic (CAFM) phase appears right above he', the 
spins acquire a uniform longitudinal component and an an- 
tiferromagnetically ordered transverse component that spon- 
taneously breaks the U(l) symmetry of global spin rotations 
along the z-axis. The CAFM phase can also be described as 
a condensation of bosonic particles. The particle density, rir, 
is related to the local magnetization along the symmetry axis 
rir = 5^ + 1. Therefore, the magnetic field acts as a chemical 
potential in the bosonic description. For hz > h^ the sys- 
tem is populated by a finite density of bosons that condense 
in the single particle state with momentum Q with = it 
{y = {x^y^z}). The longitudinal magnetization (density of 
bosons) increases with field and saturates at the fully polarized 
(FP) state (5'^ = 1 V r) above the saturation field hs. The FP 
state corresponds to a bosonic Mott insulator in the language 
of Bose gases. There exists a critical value of the single-ion 
anisotropy, Dc, below which the CAFM phase extends down 
to zero field. The nature of the QPM-CAFM quantum phase 
transition changes between hz = and hz ^ 0. The transi- 
tion belongs to the BEC universality class for hz ^ 0, while 
it belongs to the 0(2) universality class for hz = 0. 

In the next section we introduce a generalized spin wave 
theory that describes the ground state and quasiparticle ex- 
citations of the quantum paramagnetic and the canted AFM 
phases. We describe two procedures - one based on the 
standard Holstein-Primakoff approach^^, and a second one in 
which a Lagrange multiplier is introduced to enforce the local 
constraint at a mean field leveP. The QMC method is intro- 



duced in Sec.|III| Sec. [IV| includes a comparison between the 
analytical and numerical (QMC) results, which shows that the 
quantitative agreement with numerical simulations is consid- 
erably improved for the Lagrange multiplier method over the 
Holstein-Primakoff approach. We note that this is true both 
for the quantum phase diagram and for the dispersion of mag- 
netic excitations even in d = 2. This remarkable accuracy in 



describing low energy dispersion indicates that the second ap- 
proach is ideally suited for extracting Hamiltonian parameters 
from fits of INS data. SecJVjis devoted to finite temperature 
results. Finally, in Sec. |Vl[we discuss the implication of our 
results for the organic quantum magnet DTN and for any other 
quantum magnet that is close to the QCP which separates the 
magnetically ordered and paramagnetic ground states. 



II. GENERALIZED SPIN WAVE APPROACH 

In this section, we give a brief outline of the generalized 
spin wave formalism that was originally applied to the de- 
scription of the quantum paramagnetic state of DTN A Since 
the local Hilbert space has dimension Di = 3, we introduce 
three Schwinger bosons (SB) with annihilation (creation) op- 
erators bmr^m G {0, 1,2}. The three different states occu- 
pied by a single boson are mapped into the eigenstates of 
for each site r \ 

blm = \0)r, blm = \l)r, 5^10) = 1-1). (2) 
The local constraint. 



(3) 



guarantees that the dimension of the local Hilbert space is pre- 
served under this mapping. The bilinear forms of these SBs 
are generators of SU(3) in the fundamental representatioiP. 
We use the SBs to extend the usual SU(2) spin wave ap- 
proach to SU(3)^^ since the local order parameter for 5 = 1 
spins has 8 components, which correspond to the 8 gener- 
ators of the SU(3) group of unitary transformations in the 
local Hilbert space of dimension 3. Three of them corre- 
spond to the local magnetization (S*^, 5*^, 5^), while the other 
five are the components of the traceless symmetric tensor, 
Qr = {S?Sr + S!;Sll)/2 - (5^^2/3, that defines the local 
spin nematic moment. In particular, the paramagnetic mean 
field ground state has a net nematic component induced by 
the single-ion anisotropy, but no net magnetization compo- 
nent. Such a state has no classical counterpart. Nevertheless, 
we can still implement a semi-classical approximation if we 
generalize the traditional spin-wave analysis from SU(2) to 
SU(3). In this approach, we can describe the quantum fluctu- 
ations around the mean field state as small (quadratic) oscilla- 
tions of an SU(3) order parameter. 

At the mean field level, any ground state that is stabilized 
for I) > is described by the product state 



(4) 



where 



and the variational parameters and (j) are determined 
by minimization of the mean field energy per site 
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eo = {M^h\^ci)/N: 
dO 



0, 



deo 



0. 



(6) 



We note that the variational parameters 6 and (j) are enough to 
parametrize the three different phases that appear in the phase 
diagram of Tin fovD > 0. The bosonic operator b^J belongs 
to a new set of SB operators that are obtained from the original 
set by a unitary transformation, Ur- 



bor 

lArhrpi br = I bir 
b2r' 



(7) 



This transformation corresponds to choosing a quantization 
axis along the direction of the order parameter, as it is done in 
the usual spin wave treatment. Since the ground state of the 
antiferromagnetic phase breaks translational symmetry mak- 
ing the two sublattices inequivalent, the corresponding canon- 
ical transformation, Ur, is different for the two sublattices, as 
it is clear from the phase factor e*^"^ that appears in Eq.([5]). 

In terms of the SBs, the spin operators St^ assume bilinear 
forms = btS^K, 



S'r = :^Xblbor-blb2r). 

= b\^bir 



(8) 



that transform as S!^ = UrSl^Ul. The spatial dependence of 
the unitary transformation Ur can be eliminated if we change 
the original basis of the Hamiltonian 1-Lh- In particular, the 
CAFM state becomes uniform if we rotate the spin reference 
frame of one of the sublattices by angle tt along the z-axis. 
Since the uniform paramagnetic ground states oH-Lh remain 
invariant under this transformation, the unitary transforma- 
tions Ur become r -independent in the new basis for all the 
different phases of Ur^ Since, and 5^'^ ^ -5^'^, 

we have that Hh ^ 

nH = J (^vS^Sl, ^Y^^DSf -h,S',) (9) 



in the new basis, where 



1 and ar 



-1. We note 



that this change of basis shifts the AFM wave vector from Q 
to and removes the factor e*^"^ from Eq. ([5]). 

The bosonic representation of the Hamiltonian in the new 
basis is 

Uh = J a^hlS^'hrhl^S^'hr' (10) 

+ ^ ^ (l - hl^hr) -h^Y ^tS'br 



where 



S^=US^U\ A = UMA^ and A 



The condensation of the bosons b^r is implemented via a 
natural extension of the Holstein-Primakoff transformatiorP21 
to the case of more than one type of boson. From the local 
constraint ^ we obtain: 



1 - h\^hir - bl^hr 



(11) 



By applying the above condition to the Hamiltonian ([T]) and 
keeping terms up to bilinear in the bosonic creation and anni- 
hilation operators, we obtain the mean field ground state en- 
ergy 

eo = dJ^a.^o^o^oo - hj'oo + ^(1 - ^o) (12) 



and the spin wave Hamiltonian 



(r,r') 
a,/3G{l,2} 



a,i3e{l,2} 



(13) 



with the Hamiltonian parameters 

V 

Ac/3 = J y^<^i.(^aQ^go - (^Qo)^^a^) 

V 

where d is the spatial dimension. In the next step, the spin- 
wave Hamiltonian ([13]) is transformed to momentum rep- 
resentation by introducing bosonic operators in momentum 
space: 



k,a,f3 



H.C.), 
(15) 



with 



+ tal3 ^ COS kj, 

V 



(16) 



The resultant Hamiltonian can then be straightforwardly diag- 
onalized by a Bogoliubov transformation to yield the single 
particle dispersion: 



k,a ^ 



(17) 
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A. QPM phase and the Fully Polarized phase 

At the mean field level, the paramagnetic state, 
\^ci{0 = 0))=l[bl\9), 



(18) 



is the lowest energy state for large enough D, as long as 
the applied magnetic field remains below a critical value he. 
Since the unitary transformation can be chosen as the iden- 
tity, Z// = 1, the quasiparticle dispersion becomes particularly 
simple in the QPM phase: 



-2j"^cos{K). (19) 



Both branches have the same dispersion at zero field, hz = 0, 
as expected from time reversal symmetry. A finite hz splits 
the branches linearly in hz without changing the dispersion. 
This is a consequence of the fact that the external field cou- 
ples to the total magnetization, = S^, which is a con- 
served quantity. Both branches have a minimum at the AFM 
wave- vector k = that determines the size of the gap. The 
dispersion is quadratic near k = except for the critical point 
(Dc = 4(iJ, hz = 0) that separates the QPM phase from the 
CAFM phase Sithz =0. The field induced QCP then belongs 
to the BEC universality class in dimension d-\-2. By expand- 
ing around = 0, we obtain: 



ujk± ^ Jk" ^D/{D - D,) + ^D{D - D,) ± hz (20) 

It is clear from this expression that the effective mass of 
the magnetic excitations vanishes for D ^ Dc'. oc 
\/D — Dc. This is indeed the expected behavior if we keep 
in mind that the dispersion must be linear at the the critical 
point {Dc = MJ, hz = 0) (z = 1 for the 0(2) QCP as we 
discussed in the introduction). 

The QPM ground state remains stable for 



D > Dc = MJ 

hz < he = ^D{D - Dc) 



(21) 



The ground state becomes fully polarized over the saturation 
field 

hs=D + Ad J, (22) 
and the mean field state, 

|V^eK^ = V2,^ = O)) = n^L|0), (23) 

r 

coincides with the exact ground sate. The energy of the sys- 
tem is proportional to the applied field as expected. The two 
branches of magnetic excitations above the saturated state are 
given by: 



^ki hz - D - 2dJ ^rjk, 
ujk2 = 2hz- 



The flat branch, cj/c2, describes the approximated spectrum of 
two-magnon bound states that appear above a critical value of 
the single-ion anisotrop>^. 

By comparing Eqs.([T9]) and ([24]), we can see that the masses 
of the gapless bosons at the two field induced QCPs, h = he 
and h = hs can be very different: 



1 
1 

m 



Ifc^o = 2J^D/{D-De), 

|fc=o = 2 J. (25) 



While 

rate, 

rrf /m 



the mass renormalization factor rrf /m 
may not be quantitatively 



Dc)/D 

the obtained 



mean field critical 



accu- 
exponent, 

(X ^j[^u — Dc)lD, is correct for d = 3 up to 
logarithmic corrections, because (ic = 3 is the upper critical 
dimension for the 0(2) QCP in dimension d + 1. For d < dc 
we have 



mVmcx A, ex {D - Dcf 



(26) 



and the mean field exponent = 1/2 is not correct for d < 3. 
It is clear then that quantum paramagnets which are close to 
the CAFM instability {D > Dc) should exhibit a very large 
asymmetry between the mass of the bosonic excitations for 
h < he and h > hs. This is indeed the case of the compound 
DTN whose thermodynamic properties exhibit a large asym- 
metry between the the two critical points at he and hs. The 
possibility of having a relatively large m* /m ratio that can 
be tuned with pressure allows for measuring dependence of 
different physical properties on the mass of the bosonic exci- 
tations. This property of certain quantum paramagnets is par- 
ticularly useful for unveiling the dominant scattering mech- 
anism for thermal conductivity, because different mecha- 
nisms usually lead to different dependences of n on the mass 
of the quasiparticles^. 

While the linear approach that we have described gives the 
correct qualitative picture 'm d = 3, it is still far from being 
quantitatively accurate in d = 3 or (i = 2, as we will see in 
the next sections. This shortcoming can be a serious prob- 
lem for comparisons against experimental data. In particular, 
the Hamiltonian parameters for quantum paramagnets are nor- 
mally extracted from fits of the quasiparticle dispersions that 
are measured with The accuracy of the obtained Hamil- 
tonian parameters depends on the accuracy of the approach 
that is used for computing the dispersions uj^y. Moreover, 
for quantum paramagnets like DTN which have low critical 
fields, he <C hs — he, the linear approach normally predicts 
AFM ordering dXhz = 0. Therefore, it is necessary to mod- 
ify the linear approach in order to obtain a quantitatively ac- 
curate description of the low field paramagnetic ground state 
and the low-energy excitations. As we shall see in the next 
sections, the modified approach that was originally applied to 
the description of DTN^ and that we describe in the rest of this 
subsection, is quantitatively accurate for d = 3 and d = 2. 

In the modified approach we replace Eq.([TT]) by 



(24) 



(27) 
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and impose the constraint ([3j at a mean field level by intro- 
ducing the Lagrange multiplier fi\ 



mr^rnr 



(28) 



m=0 



The rest of the procedure is similar to spin-wave theory, i.e., 
we only keep terms up to quadratic order in the bosonic oper- 
ators bmr ijn = 1,2) and diagonalize the resulting quadratic 
Hamiltonian via a Bogolyubov transformation. This proce- 
dure leads to the diagonal form ([17]), but with a modified 
quasiparticle dispersion. 



(29) 



relative to the expression ([19]) that was obtained from the lin- 
ear approximation. We note that the new dispersion ( [29] ) can 
be obtained from the previous one if we replace D by ji and J 
by Js^. Therefore, in the quantum paramagnetic state, the net 
effect of including a Lagrange multiplier to enforce the con- 
straint ([3]) at the mean field level is a renormalization of the 
single-ion anisotropy and exchange parameters. 

The parameters s and ji are determined self-consistently by 
the saddle point equations^^ : 



dU 



H 



djji 



0, 



dU 



H 



ds 



0. 



(30) 



By explicitly computing the left hand side of these two equa- 
tions we obtain the following expressions: 



Vk 



-y- 

_ ^ (m + s'^r]k) 



(31) 



The stability conditions ( [2T] ) for the QPM ground state are 
replaced by 



hz < hc = ^/JiJji^^^Jic). 



(32) 
(33) 



As we will see in the next sections, the quantum phase di- 
agram that is obtained from these modified conditions is in 
much better agreement with QMC simulations. The same is 
true for the modified quasiparticle dispersion ([29]). 



B. Canted Antiferromagnetic (CAFM) phase 

To describe the CAFM phase, one needs to use the general 
expression for the condensed boson with U ^ particular, 
we use the expression given by Eq. ([5]) 

h\j, = bl^ cos 6 + sin 6> cos + 4r sin sin (j). (34) 

We recall that the factor e^^^ is removed from Eq. ([5]) after 
the change of basis that led to Eq. ([9]). The other bosonic oper- 
ators are obtained by orthogonalization. The parameters and 



are determined by the minimization of the mean field energy 
[see Eq.([6])]. In the absence of any applied field, the AFM or- 
dered phase is invariant under the product of a translation by 
one lattice parameter and a time reversal transformation. This 
symmetry implies that ^ = J, i.e., the local moments have 
equal weights in the 5'^ = ±1 states. By minimizing the 
mean field energy as a function of the remaining variational 
parameter, 0, we obtain 



1 D 



sm 



2 WdJ 



(35) 



The dispersion relation consists of two non-degenerate 
branches that, in the low energy limit (k ^ 0), are given by 



^yj{D,^D)k. 



(36) 



Unfortunately, the modified approach based on the inclu- 
sion of a Lagrange multiplier that we introduced in the previ- 
ous subsection does not work well inside the ordered phase. 
Both branches become gapped inside the ordered phase, i.e., 
the approach misses the Goldstone mode associated with the 
spontaneous breaking of the U(l) symmetry of global spin ro- 
tations along the z-axis. 

As we explained above, the magnetic field induced quan- 
tum phase transition from the QPM to the CAFM phase is 
qualitatively different from the transition between the same 
two phases that is induced by a change of D at hz = 0. 
Eq.([T9]) shows that the effect of increasing hz from zero at 
a fixed D > Dc is to reduce the gap, = uJk=o- = 
\JD^ — AdJD — hz, linearly in hz. The dispersion does not 
change because hz couples to rriz = Xlr ^r/^ that is a con- 
served quantity (m^ = 1 for the spin excitations that have dis- 
persion uJh-)- Therefore, the quasiparticle dispersion remains 
quadratic at the field induced QCP h = he = VD'^ - 4dJD, 
i.e., the dynamical exponent is z=2. The field induced QCP 
then belongs to the BEC universality class in dimension d-\-2. 
On the other hand, if the single-ion anisotropy is continuously 
decreased at zero applied field, the two branches remain de- 
generate and the gap vanishes at D = Dc (hz = 0). The low- 
energy dispersion becomes linear at the QPM-CAFM phase 
boundary, ujk± ^ y/2DJk for small k. As it is clear from 
Eq. ( [36j ), the degeneracy between the two branches dXhz = ^ 
is lifted inside the CAFM phase - one of the branches, ujk2, 
remains gapless with a linear dispersion at low energy (corre- 
sponding to the Goldstone mode of the ordered CAFM state ) 
whereas the other mode develops a gap to the lowest excita- 
tion. 

In the following sections, we shall use large scale quantum 
Monte Carlo simulations of the Hamiltonian ([TJ to demon- 
strate that the introduction of a Lagrange multiplier signif- 
icantly improves the quantitative description of the QPM 
phase, and that the linear approximation gives a qualitatively 
correct description of the quantum phase transitions in d = 3. 
As expected, md = 2, the only deviation from mean field be- 
havior occurs at the 0(2) QCP, D = Dc and hz = 0, because 
the effective dimension, P = d + 1, is lower than four. 
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IV ZERO-TEMPERATURE RESULTS 



III. QUANTUM MONTE CARLO METHOD 

We have used two different QMC methods, the standard 
stochastic series expansion (SSE) with loop update^^ and 
a modified version developed in Ref . 40 , to study the ground 
state and finite temperature properties of the Hamiltonian ([T]). 
Since both methods are unbiased and exact within the statisti- 
cal error, we refer to them as QMC collectively in this paper. 
On the dense parameter grids (temperature for thermal tran- 
sitions and magnetic field or single-ion anisotropy for ground 
state transitions) needed to study the critical region in detail, 
the statistics of the QMC results can be significantly improved 
by the use of a parallel tempering scheme^ ^ The implemen- 
tation of tempering schemes in the cont ext of the SSE method 
has been discussed in detail previousl}@^^. Ordinarily, the 
SSE would suffer from the negative sign problem for the AFM 
Heisenberg interaction. However, the sublattice rotation dis- 
cussed in section II maps the XY part of the Heisenberg in- 
teraction into a ferromagnetic exchange term, thus alleviating 
the sign problem. This transformation maps the AFM order- 
ing vector to Q = in the new basis. 

We compute the spin stiffness ps — defin ed as the response 
to a twist in the boundary conditionJ^HSI. The transition to 
CAFM is efficiently investigated by studying the scaling prop- 
erties of the spin stiffness ps. For simulations that sample 
multiple winding number sectors, the stiffness can be related 
to the fluctuations of the winding number in the updateJ^^E^EH 
and can be estimated readily with great accuracy. For the 
isotropic systems that are primarily considered in the present 
study, the estimates of the stiffness along all the axes are equal 
within statistical fluctuations. 

Along with the spin stiffness, we calculate the square of 
the order parameters characterizing the different ground states 
as well as standard thermodynamic observables such as en- 
ergy and magnetization. The transverse component of the spin 
static structure function, 

measures the off-diagonal long-range ordering in the XY 
plane. Its value at the AFM ordering wave vector, Sq~ , 
quantifies the XY AFM order. In the bosonic language, it 
is the condensate fraction of the BEC. We also compute the 
mean value of the zz-component of the nematic tensor com- 
ponent, Q^^ = {{S^Y ~ §)' that is induced by the single-ion 
anisotropy term. 



IV. ZERO-TEMPERATURE RESULTS 

A. Finite-size scaling for quantum criticality 

The continuous phase transition from the QPM phase to the 
CAFM phase is marked by the closing of the spin gap. To 
determine the transition point, we use the finite- size scaling 
properties of the spin stiffness ps . The finite-size scaling anal- 



ysis at the critical point predicts that 

below the upper critical dimension, i.e., d + z < 4, where L 
is the linear dimension of the system, z is the dynamic critical 
exponent, and Yp^ is the scaling function. 2: = 1 for QPTs 
belonging to the 0(2) universality class and z = 2 for BEC 
QCPs. Since the effective dimension of the BEC-QCP md = 
3,V = 3 + 2, is above the upper critical dimension = 4, 
we need to apply a modified finite- size scaling^^ 

The scale invariance at the critical point provides a powerful 
and widely used tool to simultaneously determine the position 
of the critical point and verify the value of z. On a plot of 
as a function of the driving parame- 
ters, D or hz, the curves for different system sizes will cross 
at the critical point provided the correct value of z is used. 

Figs. [T] shows the scaling of the stiffness close to the crit- 
ical point for the QPM-CAFM transition at hz =0 driven 
by varying the single-ion anisotropy D. From field theoretic 
arguments, the transition is expected to belong to the 0(2) 
universality class for which z = 1. Indeed, the curves were 
found to exhibit a unique crossing point only for z = 1. For 
a square lattice (top panel), we obtain a critical Dc = 5.63, in 
agreement with previous results^ ^ , whereas the transition oc- 
curs at I^c = 10.02 on a cubic lattice (bottom panel). Further 
confirmation of the 0(2) universality class of the transition is 
shown in the inset panels where on a plot of psL^^^~'^ vs. 
{D — Dc)L^^^, the data for different system sizes collapse 
onto a single curve with our estimated Dc and known critical 
exponents for the 0(2) universality class in d + 1 dimensions. 

Fig. [2] shows the modified finite- size scaling plots of the 
QPM to CAFM transition for D > Dc SiS the field is var- 
ied. The transition is expected to belong to the BEC univer- 
sality class and scale invariance for the stiffness at the critical 
point is found for z = 2 in accordance with field theoretic pre- 
dictions. Thus the analysis of the stiffness data at the quantum 
critical points show that the QPM - CAFM transition belongs 
to the 0(2) universality class for hz = 0, but changes to BEC 
universality class for hz ^ 0. 



B. Quasiparticle dispersion in the QPM phase 

The phase boundary between QPM and CAFM phases is 
also determined by the value of the single magnon excitation 
gap As. Since the Zeeman term commutes with the rest of the 
Hamiltonian, the spin gap of the QPM phase changes linearly 
in the magnetic field and vanishes at the critical field he = 
As{hz = 0). The quasiparticle dispersion and the gap A5 can 
be extracted from the QMC results by analysing the imaginary 
time Green's function 

GTir) = ^J2(^rir)Sm)e"'-. (38) 
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FIG. 1: (Color online) Finite-size scaling plots of spin stiffness ps- 
The four system sizes of the square lattices (upper panel) L x L are 
8x8 (red), 10 x 10 (blue), 12 x 12 (black) and 18 x 18 (purple). 
The five system sizes of the cubic lattices (lower panel) L x L x L 
are 4 X 4 X 4 (red), 6x6x6 (blue), 8x8x8 (black), 10 x 10 x 10 
(purple) and 12x12x12 (green). The temperatures are taken to 
be T = L/4 in the square lattice and L/2 in the cubic lattice. The 
boundary conditions are periodic. 



The quasiparticle dispersion is computed by fitting the QMC 
data of G^{t) with the function 



/(r) = A e 



-u;(/3-r) 



(39) 



where A and uo are fitting parameters. In particular, the pa- 
rameter uo corresponds to the magnetic excitation energy for 
each momentum k. Figure|3]shows that the fit is nearly perfect 
for the G%^{r) curve that is obtained in the QPM phase. The 
estimated phase boundary is he = 4.2726(3) for D = 12, 
d = 3 and L = 12. This estimation is fully consistent 
with the modified finite-size scaling analysis. (See Fig. [2]) 
Since finite size effects are very small deep inside the QPM 
state (far from critical point), the field induced phase bound- 
ary can be estimated very precisely with L = 12. Fig. |4] 
shows the comparison between the quasiparticle dispersions 
obtained from the QMC results and the analytical expressions 



D=12.0, h =4.273 




%.00 -30.00 -20.00 -10.00 0.00 10.00 20.00 



FIG. 2: (Color online) Determination of the critical field through 
finite size scaling with z = 2 that confirms the BEC universality 
class of the field induced quantum critical points. 




FIG. 3: (Color online) Imaginary time Green's function computed 
with QMC for D = 12, and hz = 0. The linear size of the finite 
cubic lattice is L = 12 and the boundary conditions are periodic. 
The solid fitting lines correspond to the function defined in Eq. |39l. 



(T9\ and ( |29| ) that we derived in the previous section using the 
Holstein-Primakoff (HP) and the Lagrange multiplier (LM) 
approaches. The quantitative agreement with the numerical 
result is much better for the LM approach that reproduces not 
only the value of the spin gap and the overall dispersion in- 
side the QPM phase, but also the spin velocity at the 0(2) 
QCF D = Dc{h^ = 0). 



C. Quantum phase diagram 

The quantum phase diagrams obtained with different meth- 
ods: linear HP approximation, the LM approach and QMC 
simulations, are shown in Figs. |5] As it is expected from the 
comparisons between the quasiparticle dispersions obtained 
with the different methods in the QPM phase (see Fig.|4]), the 
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FIG. 4: (Color online) Dispersions of the single magnon excitation 
(SL) D = 8 in 2D, (b) D = Dc in 2D and (c) i:) = 12 in 3D. In 2D, 
Dc = 8, 5.71 and 5.625 for the linear HP, LM and QMC approaches, 
respectively. 




FIG. 5: (Color online) Quantum phase diagram of Hh in (a) = 2 
and (b) d = 3. The solid line, dashed line and points between QPM 
and CAFM are the results obtained from the LM, HP and QMC ap- 
proaches, respectively. For the QMC approach we use the modified 
finite- size scaling that is described in the text as well as the gap that 
is obtained from the quasiparticle dispersion to determine the QPM- 
CAFM phase boundary. 



LM method produces a much better quantitative agreement 
with the QMC results than the linear HP approximation. 

Fig.[6]shows the evolution of some observables that charac- 
terize the ground state phases as the applied field is varied 
for three representative values of the single-ion anisotropy. 
For D > Dc, the ground state evolves from a QPM phase 
at low fields (h^ < he) to a. CAFM phase at intermediate 
fields (he < hz < /is) to a fully polarized phase at large 
fields. The uniform magnetization, m^, increase monoton- 
ically with the applied field. The zz-nematic order param- 
eter, Q^^, also increases monotonically but from a negative 
to a positive value. Right above h = h^ the magnetization 
rriz increases with finite slope, but this slope vanishes at the 
0(2) QCP where hc{Dc) = 0. This result is consistent with 
the mean field theory described in the previous section which 
predicts that ex {hz — hc{D)) for finite hc{D) and small 
enough hz — hc{D), while cx h^ for he = and small 
enough hz. These results are obtained by solving Eqs.([6]) near 
the 0(2) QCP {D = Dc.hz = 0). 

The stiffness and transverse structure factor decrease mono- 
tonically with increasing hz for D ^ Dc. However, it is clear 



that the field dependence must be non-monotonic for D > D^ 
because a finite critical field is required to induce the transition 
from the QPM to the ordered XY phase. When the system is 
in the QPM phase, a critical field hc{D) is required to induce 
a finite amplitude of the XY order parameter, i.e., the mean 
field state of each spin becomes a linear combination of the 
states 10)-^ and |1)^ for > he. There is an optimal value 
of the magnetic field, hm{D), for which the weight of these 
two states is roughly the same, leading to maxima of the order 
parameter (XY component of the local moment) and the spin 
stiffness, as it is shown in Fig. [6] Finally, ps and S^-{Q) 
vanish again at sufficiently strong applied field, hz > hs{D), 
because the ground state evolves to the fully polarized phase 
with rriz = I, and Q^^ = 1/3. The exact boundary between 
the CAFM and the FP phases is given by Eq.([22|). A simple 
continuity argument shows that the non-monotonic field de- 
pendence of Ps and (Q) should persists for D < Dc sls 
it is clear from Fig. [6] The ordering temperature should also 
exhibit a similar non-monotonic field dependence, as we will 
see in the next section. This observation can be used to detect 
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FIG. 6: (Color online) The evolution of various characteristic ob- 
servables with external magnetic field at three representative values 
of D as the ground state goes through the field driven quantum phase 
transitions discussed in the text. The data is for a finite cubic lattice 
of dimension 16 x 16 x 16. 



quantum magnets that exhibit magnetic ordering at /i^ = , 
but are near the QCP, i.e., close to becoming quantum param- 
agnets. 



V. FINITE-TEMPERATURE RESULTS 




FIG. 7: (Color online) The critical temperatures of the thermal phase 
transition into different ground states shown in Fig.|5jb). 

For three-dimensional systems, the CAFM phase survives 
up to a finite temperature Tc{D^ hz) above which the system 
becomes a paramagnet via a second order classical phase tran- 
sition that belongs to the 0(2) universality class in dimension 
d. The second order transition is replaced by a Berezinskii- 
Kosterlitz-Thouless phase transition at T = Tbkt when the 



system is two-dimensional. In this case, only quasi long range 
ordering survives at finite temperatures T < Tbkt- Fig-[2] 
shows the field dependence of the critical temperature, Tc, for 
some representative values of D. Tc is determined by exploit- 
ing the scale invariance of the stiffness at the critical point 
with the finite-size scaling 



p,(L,T) ~ L^-%,((T-Te)L^/'^). 

The thermal transition out of the CAFM phase is driven by 
phase fluctuations of the order parameter and belongs to the 
(i = 3 0(2) universality class {y 0.67). At small values 
of D, the system is dominated by the Heisenberg AFM in- 
teraction and Tc{hz) decreases monotonically as a function 
of increasing hz to Tc{hs) =0 at the QMP-FP boundary. 
As D increases, the spins acquire a significant 5*^ = (ne- 
matic) component and the resultant decrease in the local mag- 
netization leads to a suppression of the critical temperature. 
As we explained in the previous section, the applied field in- 
creases the magnitude of the local moments for D < Dc and 
this effect leads to an accompanying increase in Tc{h). At 
higher values of the applied field, the spins acquire an increas- 
ing (ferromagnetic) component along the field direction while 
the AFM-ordered component decreases beyond the optimal 
field hm{D). Consequently, the critical temperature starts de- 
creasing monotonically to Tc{hs) = for h > hm{D). For 
D > Dc, the system is in a QPM ground state at low fields - 
with the local spins being predominantly in the Sz = state 
- and Tc = 0. A sufficiently strong external field induces a 
transition to the CAFM phase with Tc cx {hz - hcY^^ for 
small enough hz — he. The transition temperature increases 
initially as the magnitude of the local moments increase and 
eventually decreases as the moments acquire a dominant fer- 
romagnetic component parallel to the applied field - going to 
Tc = OsLthz = hs. 



VI. SUMMARY 

In summary, we have investigated the quantum phase di- 
agram and the nature of the quantum phase transitions in 
the 5 = 1 Heisenberg model with easy-plane single-ion 
anisotropy and an external magnetic field. By using a gen- 
eralized spin wave approach, we showed that the low energy 
quasiparticle dispersion is qualitatively different at the phase 
boundary depending on the presence or absence of an exter- 
nal field. This difference is reflected in the universality class 
of the underlying QCP and has direct consequences on the low 
temperature behavior. The nature of the QPM-CAFM transi- 
tion in the presence and absence of an external field is directly 
confirmed by using large scale QMC simulations and finite 
size scaling. 

We have used two different analytical approaches to de- 
scribe the QPM. By comparing the results of both approaches 
against our QMC results, we have found important quanti- 
tative differences in the region near the 0(2) QCP that sig- 
nals the transition to the CAFM phase. By "quantitative dif- 
ferences" we are not referring to the already known criti- 
cal behaviors predicted by both approaches, but to the phase 
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boundary Dc{hz) and the dispersion of the low-energy quasi- 
particle excitations. To make a clear distinction between these 
two different aspects of the problem, we will discuss the crit- 
ical behavior in the first place. It is clear that both analytical 
treatments reproduce the correct critical behavior for d = 3 
up to logarithmic corrections, because dc > 3 for the QCPs 
[0(2) and BEC] that appear in the quantum phase diagram of 
Hh- The situation is different for d = 2 because the upper 
critical dimension of the 0(2) QCP is dc = 3. We note that 
the approach based on the inclusion of the Lagrange multi- 
plier and the saddle point approximation ( [3Q| , becomes exact 
in the large A/" ^ oo limit (A/" is the number of components of 
the order parameter of the broken symmetry state, i.e.. A/" = 2 
for the case under consideration)^^. Since = l/((i — 1) for 
A/" cx), the LM approach leads to a spin gap that closes 
linearly in {D — Dc) for d = 2 (see Fig. |5^). In contrast, 
the HP approach produces the expected mean field exponent 
v = 1/2. Naturally, neither of these approaches can repro- 
duce the correct value of the exponent u [v 0.67 for the 
0(2) QCP in dimension P = 2 + 1] because 2 < dc. However, 
the LM approach can be systematically improved by including 
higher order corrections in l/A/". The qualitative agreement 
for d = 3 is not surprising because the effective dimension of 
the QCPs that appear in the quantum phase diagram of Hh is 
equal or higher than the upper critical dimension. 

Since the limitations of the LM and HP approaches for de- 
scribing the critical behavior of the 0(2) QCP are already 
known, we have focused on the overall quantitative agree- 
ment for the phase boundary Dc{hz) and the dispersion of 
the low-energy quasi-particle excitations in comparison with 
the numerical results. The very good agreement between the 
LM and QMC results is rather surprising if we consider that 
it holds true even for (i = 2 (see Fig. |4] and |5^). Indeed, a 
similar treatment has been successfully applied to the quasi- 
one-dimensional organic quantum magnet known as DTN ^ 



In this compound, the S = 1 moments are provided by Ni^+ 
ions which are arranged in a tetragonal lattice. The magnetic 
properties are well described by the Hamiltonian ([T]) with pa- 
rameters D S.9K, Jc 2.2K and Ja = Jb = O.ISK, 
where Ja denotes the strength of the Heisenberg exchange 
interaction along the different crystal axes. Once again, the 
introduction of a Lagrange multiplier to enforce the constraint 
^ leads to a critical field value of 2T, that is in very good 
agreement with the result of QMC simulations and with the 
experiments^ ^. In contrast, the linear HP approach incorrectly 
predicts that this compound should be magnetically ordered 
in absence of the applied magnetic field. We note that the 
phase boundary obtained with the LM approach fovd = 2 
(see Fig. |5^) remains quantitatively more accurate near the 
0(2) QCP even when the next (second) order corrections in 
1/S are included in the HP approach^^. Our results then indi- 
cate that introducing a Lagrange multiplier for describing the 
low-energy physics of quantum paramagnets improves con- 
siderably the estimation of the spin gap and the quasiparti- 
cle dispersion. This improvement is particularly important for 
quantum paramagnets that have a small spin gap and conse- 
quently are close to the QCP that signals the onset of magnetic 
ordering. Since the Hamiltonian parameters are typically ex- 
tracted from fits of the quasiparticle dispersion measured with 
INS, it is crucial to have a reliable approach for computing 



such dispersion. The QMC method described in Sec. |IVB| can 
only be applied to Hamiltonians that are free of the sign prob- 
lem. However, the analytical approach described in Sec.[II|is 
always applicable. 

The numerical results were obtained in part using the com- 
putational resources of the National Energy Research Scien- 
tific Computing Center, which is supported by the Office of 
Science of the U.S. Department of Energy under Contract No. 
DE-AC02-05CH11231. 
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